Setting

path_proj = here::here()
path_source = file.path(path_proj, "source")

source(file.path(path_source, "simulation", "simulations_functions_final.R"))
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:reshape2':
## 
##     smiths
source(file.path(path_source, "functions", "plot_function.R"))
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
source(file.path(path_source, "functions", "fit_function.R"))
source(file.path(path_source, "functions", "table_function.R"))

# place for draws
# mac
posterior_draws_path = file.path(Sys.getenv("HOME"), "Desktop", "draws", "testEach")

#Windows
#posterior_draws_path = file.path(Sys.getenv("USERPROFILE"), "Desktop", "draws", "testEach")

#data path
data_save_path = file.path(path_proj, "data", "fitted_model", "simulation", "3. b_rw")
#models
q_constant <- file.path(path_proj, "source", "models", 
                      "q-constant.stan")
b_constant <- file.path(path_proj, "source", "models", 
                     "b-constant.stan")
b_rw <- file.path(path_proj, "source", "models", 
                     "b-rw1.stan")
b_ou <-  file.path(path_proj, "source", "models", 
                      "b-ou.stan")

compiled_models <- list(
  q_constant = cmdstan_model(q_constant),
  b_constant = cmdstan_model(b_constant),
  b_rw = cmdstan_model(b_rw),
  b_ou = cmdstan_model(b_ou)
)

models_to_use <- c("q_constant", "b_constant", "b_rw", "b_ou")
###### setting #####
# data
alpha_increase_seq_1 <- seq(10, 750, by = 30)
alpha_decrease_seq_1 <- seq(750, 10, by = -30)
alpha_lamb =  c( rep(10,5), alpha_increase_seq_1 + rnorm(alpha_increase_seq_1,10,10), 
                 alpha_decrease_seq_1 + rnorm(alpha_decrease_seq_1,10,10),
                 rep(10,5))
beta_lamb = 0.5
T = 60
# reprot delay
D <- 15;

# Time period for checking
D_check <- 5

first_date <- as.Date("2024-01-01")

scoreRange <- c(first_date+days(9), first_date+days(19), first_date+days(29),
                first_date+days(39), first_date+days(49))

Fully reported case

Simulation

params_b_rw_FR <- list(
  data = list(
    alpha_lamb = alpha_lamb,
    beta_lamb  = beta_lamb,
    T       = T,
    date_start = as.Date("2024-01-01"),
    D = D
  ),
  q_model = list(
    method        = "b_rw",
    method_params = list(mu_logb = log(0.7), sigma_logb = 0.1, mu_logitphi = 1, sigma_logitphi = 0.1)
  )
)

b_rw_FR <- simulateData(params_b_rw_FR)
par(mfrow = c(2, 1))
plot(b_rw_FR$b, pch = 19, type = "b")
plot(b_rw_FR$phi, pch = 19, type = "b")

par(mfrow = c(1, 1))
matplot(t(b_rw_FR$q), type = "l", lty = 1, ylim = c(0, 1))

Exploratory analysis

# exploritary analysis
page_num <- ceiling(nrow(b_rw_FR$case_reported_cumulated)/16)
exp_plot_b_rw <- fit_exp_plot(b_rw_FR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# print(exp_plot_b_rw)
# 
# exp_plot_b_rw$coefficients
exp_b_data_b_rw<- data.frame( date = as.Date(rownames(b_rw_FR$case_reported_cumulated)),
                          b = exp_plot_b_rw$coefficients$b)

exp_b_plot_b_rw <- ggplot(exp_b_data_b_rw, aes(x = date, y = b)) +
  geom_point(color = "black", size = 1.5) +       
  geom_smooth(method = "loess", se = TRUE,        
              color = "blue", fill = "grey", alpha = 0.5) +
  theme_minimal() +
  labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")

print(exp_b_plot_b_rw)
## `geom_smooth()` using formula = 'y ~ x'

Model fitting

# out_b_rw_FR <-
#   nowcasting_moving_window(b_rw_FR$case_reported_cumulated, scoreRange = scoreRange,
#                           case_true = b_rw_FR$case_true,
#                           start_date = as.Date("2024-01-01"),
#                           D = D,
#                           methods =models_to_use,
#                           compiled_models = compiled_models,
#                           iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
#                           num_chains = 3, suppress_output = T,
#                           posterior_draws_path = file.path(posterior_draws_path, "b_rw")
#                           )
# 
# save(out_b_rw_FR, file = file.path(data_save_path, "FR_b_rw.RData"))
load( file.path(data_save_path, "FR_b_rw.RData"))

results_b_rw_FR <- nowcasts_table(out_b_rw_FR, D = D, report_unit = "day", 
                          methods = models_to_use)

plots_b_rw_FR <- nowcasts_plot(results_b_rw_FR, D = D, report_unit = "day", methods = models_to_use)
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
for (i in 1:length(results_b_rw_FR)) {
  print(calculate_metrics(results_b_rw_FR[[i]],
                          methods = models_to_use))
}
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 77.76 62.00 53.53 60.54         120.61           0.6 q_constant
## 2 15.48  8.05  8.07  6.66          57.31           1.0 b_constant
## 3 16.19  8.09  8.20  6.71          67.21           1.0       b_rw
## 4 16.11  8.18  8.19  6.75          61.61           1.0       b_ou
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 67.30 12.43 46.87 11.57         107.16          0.75 q_constant
## 2 31.90  4.50 15.75  3.07          80.65          1.00 b_constant
## 3 27.66  3.88 12.12  2.40          95.06          1.00       b_rw
## 4 32.36  4.49 15.30  2.93          90.51          1.00       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 54.52  4.65 34.19 3.75         109.34          0.87 q_constant
## 2 30.63  2.88 15.09 2.07         103.67          0.97 b_constant
## 3 22.12  1.91  8.14 1.12         124.34          1.00       b_rw
## 4 22.88  2.04  9.28 1.27         124.95          1.00       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 24.40  2.43 15.79 1.82         114.33             1 q_constant
## 2 17.62  2.33  9.94 1.59         112.18             1 b_constant
## 3 19.45  2.14  7.97 1.10         126.29             1       b_rw
## 4 13.36  1.71  6.42 1.02         125.13             1       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 33.60  6.90 15.32 2.61         110.07          0.98 q_constant
## 2 30.66  6.63 11.12 2.32         108.13          0.98 b_constant
## 3 15.13  3.31  5.80 1.24         112.65          1.00       b_rw
## 4 19.68  4.35  6.14 1.40         113.17          1.00       b_ou
for (i in 1:length(results_b_rw_FR)) {
  print(calculate_metrics(data.table::last(results_b_rw_FR[[i]],D_check),
                          methods = models_to_use))
}
##     RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 109.38 67.98 96.61 66.08         211.23           0.2 q_constant
## 2  21.86  9.17 15.36  8.42          97.01           1.0 b_constant
## 3  22.86  9.31 15.58  8.40         116.41           1.0       b_rw
## 4  22.75  9.33 15.60  8.43         105.81           1.0       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 117.55 15.88 109.96 14.80         226.63           0.2 q_constant
## 2  61.62  7.75  49.53  6.40         166.02           1.0 b_constant
## 3  53.39  6.65  37.18  4.76         219.22           1.0       b_rw
## 4  62.92  7.89  49.07  6.33         202.02           1.0       b_ou
##     RMSE RMSPE    MAE MAPE Interval_Width Coverage_Rate     Method
## 1 110.77  8.05 102.23 7.52         231.61           0.4 q_constant
## 2  68.07  4.85  58.24 4.25         220.60           0.8 b_constant
## 3  52.28  3.67  34.30 2.43         337.63           1.0       b_rw
## 4  54.72  3.83  42.74 3.04         339.24           1.0       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 46.65  4.58 35.68 3.39         198.42             1 q_constant
## 2 38.92  3.77 28.06 2.64         192.42             1 b_constant
## 3 53.08  5.24 39.78 3.83         293.43             1       b_rw
## 4 35.47  3.46 28.97 2.77         289.23             1       b_ou
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 97.61 21.08 62.40 13.01         140.01           0.8 q_constant
## 2 93.36 20.19 56.80 11.93         137.80           0.8 b_constant
## 3 45.81  9.84 31.14  6.40         169.40           1.0       b_rw
## 4 61.00 13.19 38.42  7.95         175.02           1.0       b_ou
print(plots_b_rw_FR)
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

Parameter Checking

# try the third one, "2024-01-30"
T_test = T * 3/6

# q_constant
varnames_b_rw <- out_b_rw_FR$b_rw[[3]]$summary()$variable


mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws("sigma_log_b"), prob_outer = 0.95)

mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws("sigma_logit_phi"), prob_outer = 0.95)

param_true = tibble(
    parameter = grep("^b\\[.+\\]$", varnames_b_rw, value = TRUE),
    x = b_rw_FR$b[1:T_test]
)
mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws("b"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).

param_true = tibble(
    parameter = grep("^phi\\[.+\\]$", varnames_b_rw, value = TRUE),
    x = b_rw_FR$phi[1:T_test]
)

mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws("phi"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).

param_true = tibble(
    parameter = grep("^lambda\\[.+\\]$", varnames_b_rw, value = TRUE),
    x = b_rw_FR$lambda[1:T_test]
)
mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws("lambda"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^q\\[10,.+\\]$", varnames_b_rw, value = TRUE),
    x = b_rw_FR$q[10,]
)
mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws(grep("^q\\[10,.+\\]$", varnames_b_rw, value = TRUE)), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

param_true = tibble(
    parameter = grep("^N\\[.+\\]$", varnames_b_rw, value = TRUE),
    x = b_rw_FR$case_true[1:T_test, 1]
)
mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws("N"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

Simulation

params_b_rw_NFR <- list(
  data = list(
    alpha_lamb = alpha_lamb,
    beta_lamb  = beta_lamb,
    T       = T,
    date_start = as.Date("2024-01-01"),
    D = D
  ),
  q_model = list(
    method        = "b_rw",
    method_params = list(mu_logb = log(0.1), sigma_logb = 0.1, mu_logitphi = 1.5,
            sigma_logitphi = 0.1)
  )
)

b_rw_NFR <- simulateData(params_b_rw_NFR)
par(mfrow = c(2, 1))
plot(b_rw_NFR$b, pch = 19, type = "b")
plot(b_rw_NFR$phi, pch = 19, type = "b")

par(mfrow = c(1, 1))
matplot(t(b_rw_NFR$q), type = "l", lty = 1, ylim = c(0, 1))

b_rw_NFR$case_true - b_rw_NFR$case_reported_cumulated[,16]
##            [,1]
## 2024-01-01    1
## 2024-01-02    1
## 2024-01-03    3
## 2024-01-04    1
## 2024-01-05    3
## 2024-01-06    1
## 2024-01-07   13
## 2024-01-08   11
## 2024-01-09   17
## 2024-01-10   28
## 2024-01-11   36
## 2024-01-12   41
## 2024-01-13   41
## 2024-01-14   43
## 2024-01-15   41
## 2024-01-16   50
## 2024-01-17   41
## 2024-01-18   57
## 2024-01-19   57
## 2024-01-20   59
## 2024-01-21   69
## 2024-01-22   79
## 2024-01-23  100
## 2024-01-24  145
## 2024-01-25  131
## 2024-01-26  194
## 2024-01-27  217
## 2024-01-28  269
## 2024-01-29  250
## 2024-01-30  165
## 2024-01-31  235
## 2024-02-01  192
## 2024-02-02  222
## 2024-02-03  357
## 2024-02-04  377
## 2024-02-05  339
## 2024-02-06  380
## 2024-02-07  304
## 2024-02-08  348
## 2024-02-09  378
## 2024-02-10  281
## 2024-02-11  281
## 2024-02-12  260
## 2024-02-13  268
## 2024-02-14  263
## 2024-02-15  184
## 2024-02-16  168
## 2024-02-17  143
## 2024-02-18  133
## 2024-02-19  116
## 2024-02-20  106
## 2024-02-21   95
## 2024-02-22   65
## 2024-02-23   61
## 2024-02-24   25
## 2024-02-25    6
## 2024-02-26    1
## 2024-02-27    4
## 2024-02-28   10
## 2024-02-29    4

Exploratory analysis

# exploritary analysis
page_num <- ceiling(nrow(b_rw_NFR$case_reported_cumulated)/16)
exp_plot_b_rw <- fit_exp_plot(b_rw_NFR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
# print(exp_plot_b_rw)
# 
# exp_plot_b_rw$coefficients
exp_b_data_b_rw<- data.frame( date = as.Date(rownames(b_rw_NFR$case_reported_cumulated)),
                          b = exp_plot_b_rw$coefficients$b)

exp_b_plot_b_rw <- ggplot(exp_b_data_b_rw, aes(x = date, y = b)) +
  geom_point(color = "black", size = 1.5) +       
  geom_smooth(method = "loess", se = TRUE,        
              color = "blue", fill = "grey", alpha = 0.5) +
  theme_minimal() +
  labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")

print(exp_b_plot_b_rw)
## `geom_smooth()` using formula = 'y ~ x'

Model fitting

# out_b_rw_NFR <-
#   nowcasting_moving_window(b_rw_NFR$case_reported_cumulated, scoreRange = scoreRange,
#                           case_true = b_rw_NFR$case_true,
#                           start_date = as.Date("2024-01-01"),
#                           D = D,
#                           methods =models_to_use,
#                           compiled_models = compiled_models,
#                           iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
#                           num_chains = 3, suppress_output = T,
#                           posterior_draws_path = file.path(posterior_draws_path, "b_rw")
#                           )
# 
# save(out_b_rw_NFR, file = file.path(data_save_path, "NFR_b_rw.RData"))
load( file.path(data_save_path, "NFR_b_rw.RData"))

results_b_rw_NFR <- nowcasts_table(out_b_rw_NFR, D = D, report_unit = "day", 
                          methods = models_to_use)

plots_b_rw_NFR <- nowcasts_plot(results_b_rw_NFR, D = D, report_unit = "day", methods = models_to_use)

for (i in 1:length(results_b_rw_NFR)) {
  print(calculate_metrics(results_b_rw_NFR[[i]],
                          methods = models_to_use))
}
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 25.18 13.66 13.34 11.68          83.11           1.0 q_constant
## 2 81.13 50.91 53.68 50.29          34.20           0.1 b_constant
## 3 81.50 50.76 53.87 50.10          35.71           0.1       b_rw
## 4 80.31 50.53 53.15 49.92          35.60           0.1       b_ou
##     RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 119.63 19.24 75.22 16.97          92.81          0.60 q_constant
## 2 129.31 21.16 83.42 19.22          82.61          0.55 b_constant
## 3  66.77 11.81 38.68  9.65         175.16          1.00       b_rw
## 4  76.49 13.18 40.90 10.63         170.07          1.00       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 183.66 17.40 121.32 15.64         101.94          0.40 q_constant
## 2 122.54 10.38  65.11  7.48         112.74          0.83 b_constant
## 3 176.11 14.76  85.15 10.05         177.45          0.87       b_rw
## 4 187.74 15.32  89.54 10.13         160.38          0.80       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 249.41 21.24 180.41 18.36         103.21          0.30 q_constant
## 2 154.88 13.37 106.39 10.93         116.46          0.55 b_constant
## 3  85.01  9.56  49.94  6.37         231.75          0.98       b_rw
## 4 126.04 12.41  73.58  7.98         209.96          0.88       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 212.52 20.39 159.93 17.66          99.91          0.30 q_constant
## 2 108.60 12.29  87.94 11.30         112.79          0.36 b_constant
## 3  41.08  7.25  28.40  4.96         195.04          1.00       b_rw
## 4  54.01  8.99  35.38  5.77         188.19          0.98       b_ou
for (i in 1:length(results_b_rw_NFR)) {
  print(calculate_metrics(data.table::last(results_b_rw_NFR[[i]],D_check),
                          methods = models_to_use))
}
##     RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1  35.55 15.73 24.95 13.57         143.22             1 q_constant
## 2 114.37 54.67 98.70 54.24          55.01             0 b_constant
## 3 114.90 54.79 99.12 54.35          58.01             0       b_rw
## 4 113.21 54.24 97.68 53.79          57.80             0       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 227.40 28.63 219.14 28.19         183.62             0 q_constant
## 2 244.33 30.84 236.68 30.49         164.01             0 b_constant
## 3 129.16 16.04 120.63 15.40         443.04             1       b_rw
## 4 149.88 18.07 130.84 16.39         423.83             1       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 395.59 28.62 392.93 28.30         209.42           0.0 q_constant
## 2 284.81 20.75 276.31 20.01         238.01           0.2 b_constant
## 3 423.53 30.19 418.58 29.90         498.06           0.2       b_rw
## 4 452.63 32.16 446.30 31.82         423.05           0.0       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 335.88 31.72 332.02 30.88         169.41           0.0 q_constant
## 2 216.86 21.23 208.19 19.64         194.81           0.0 b_constant
## 3 189.36 19.85 168.41 16.52         575.30           0.8       b_rw
## 4 287.13 28.59 277.55 26.49         423.63           0.4       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 152.48 28.22 147.75 27.96         119.22           0.0 q_constant
## 2  81.47 14.80  75.28 14.10         138.21           0.6 b_constant
## 3  72.75 14.07  69.27 13.35         326.43           1.0       b_rw
## 4 116.70 22.80 114.97 22.30         252.42           0.8       b_ou
print(plots_b_rw_NFR)
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

Parameter Checking

# try the third one, "2024-01-30"
T_test = T * 3/6

# q_constant
varnames_b_rw <- out_b_rw_NFR$b_rw[[3]]$summary()$variable


mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws("sigma_log_b"), prob_outer = 0.95)

mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws("sigma_logit_phi"), prob_outer = 0.95)

param_true = tibble(
    parameter = grep("^b\\[.+\\]$", varnames_b_rw, value = TRUE),
    x = b_rw_NFR$b[1:T_test]
)
mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws("b"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

param_true = tibble(
    parameter = grep("^phi\\[.+\\]$", varnames_b_rw, value = TRUE),
    x = b_rw_NFR$phi[1:T_test]
)

mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws("phi"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^lambda\\[.+\\]$", varnames_b_rw, value = TRUE),
    x = b_rw_NFR$lambda[1:T_test]
)
mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws("lambda"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^q\\[10,.+\\]$", varnames_b_rw, value = TRUE),
    x = b_rw_NFR$q[10,]
)
mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws(grep("^q\\[10,.+\\]$", varnames_b_rw, value = TRUE)), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

param_true = tibble(
    parameter = grep("^N\\[.+\\]$", varnames_b_rw, value = TRUE),
    x = b_rw_NFR$case_true[1:T_test, 1]
)
mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws("N"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)